library(tidyverse)
library(SingleCellExperiment)
library(Statial)
library(spicyR)
library(plotly)
library(ClassifyR)
library(lisaClust)
library(ggsurvfit)
library(ggthemes)
library(DT)
library(htmlwidgets)
library(ggthemes)
# Set ggplot theme
theme_set(theme_classic())
Dataset 1 - Keren et al
Dataset 2 - Ali et al
The next sections will look at Keren et al.
load("../data/kerenSCE.rda")
#load("../data/aliSCE.rda")
kerenSCE
## class: SingleCellExperiment
## dim: 48 197678
## metadata(0):
## assays(1): intensities
## rownames(48): Na Si ... Ta Au
## rowData names(0):
## colnames(197678): 1 2 ... 197677 197678
## colData names(42): x y ... Censored region
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
# Removing patients with cold tumour due to small sample size of 5
kerenSCE = kerenSCE[, kerenSCE$tumour_type != "cold"]
kerenSCE$tumour_type = droplevels(kerenSCE$tumour_type)
# Removing patients 22 and 38 because they have missing survival data.
kerenSCE = kerenSCE[, !kerenSCE$imageID %in% c("22", "38")]
kerenSCE$imageID = as.numeric(kerenSCE$imageID)
# Extracting clinical information
clinicalDf = kerenSCE |>
colData() |>
data.frame() |>
select(-c(x, y, CellID, cellType, cellSize, C, tumorYN, tumorCluster, Group, immuneCluster, immuneGroup, region)) |>
unique() |>
remove_rownames()
set.seed(51773)
kerenSCE <- scater::runUMAP(kerenSCE, exprs_values = "intensities", name = "UMAP")
scater::plotReducedDim(kerenSCE, dimred = "UMAP", colour_by = "cellType")
p = kerenSCE |>
colData() |>
data.frame() |>
filter(imageID == "5") |>
ggplot(aes(x = x, y = y, col = cellType)) +
geom_point(size = 1) +
ggthemes::scale_colour_tableau( palette = "Tableau 20")
ggplotly(p)
cellProp = spicyR::getProp(
cells = kerenSCE
)
cellProp |>
round(4) |>
DT::datatable(options = list(scrollX = TRUE))
testProp = spicyR::colTest(kerenSCE,
condition = "tumour_type",
feature = "cellType")
testProp |>
DT::datatable(options = list(scrollX = TRUE))
cellProp |>
rownames_to_column("imageID") |>
mutate(imageID = as.numeric(imageID)) |>
left_join(clinicalDf) |>
ggplot(aes(x = tumour_type, y = `Keratin+Tumour`)) +
geom_boxplot() +
labs(y = "Proprotion of Keretin+Tumour cells",
x = "Tumour Type")
set.seed(51773)
kerenSCE <- lisaClust(
kerenSCE,
k = 5,
spatialCoords = c("x", "y"),
cellType = "cellType"
)
regionMap(kerenSCE,
cellType = "cellType")
hatchingPlot(
kerenSCE,
useImages = "5",
cellType = "cellType",
spatialCoords = c("x", "y"),
line.spacing = 41, # spacing of lines
nbp = 100 # smoothness of lines
)
A positive value means the two cell types are localised, a negative value means they are dispered.
interactions = spicyR::getPairwise(
cells = kerenSCE,
) |>
data.frame()
interactions |>
round(4) |>
DT::datatable( options = list(scrollX = TRUE))
The interaction between Keratin+Tumour cells and other immune cells are the most significant spatial relationship which contributes to patient survival.
The CoxPH coefficient is negative, which means when Keratin+Tumour and other immune cells are localising, the patient tends to live longer.
# Create survival object from data
survivalOutcomes = Surv(clinicalDf$`Survival_days_capped.`, clinicalDf$Censored)
# Fit CoxPh models on all cell relationships.
ClassifyR::colCoxTests(interactions, survivalOutcomes) |>
arrange(p.value) |>
round(4) |>
DT::datatable( options = list(scrollX = TRUE))
# Selecting the most significant relationship
relationship = interactions$Keratin.Tumour__other.immune
# Splitting the values by median relationship
relationship = ifelse(relationship > median(relationship), "Attraction", "Avoidance")
# Plotting Kaplan-Meier curve
survfit2(survivalOutcomes ~ relationship) |>
ggsurvfit() +
add_pvalue() +
ggtitle("Keratin+Tumour__Other immune")
# Calculate proportions of regions in each image
regionProp <- spicyR::getProp(kerenSCE,
feature = "region",
imageID = "imageID")
# Get the average expression of a marker in each cell type in each region
cellTypeRegionMeans <- Statial::getMarkerMeans(kerenSCE,
imageID = "imageID",
cellType = "cellType",
region = "region")
# Group all features in one vector
featureList = list(
"Cell type proportions" = cellProp,
"Region proportions" = regionProp,
"Marker means" = cellTypeRegionMeans,
"Spatial interactions" = interactions
)
# Training a CoxNet survival model
CV = ClassifyR::crossValidate(
measurements = featureList,
outcome = survivalOutcomes,
classifier = "CoxNet",
selectionMethod = "CoxPH",
nFolds = 5,
nFeatures = 10,
nRepeats = 20
)
ClassifyR::performancePlot(CV, metric = "C-index",
characteristicsList = list(x = "auto", fillColour = "Assay Name")) +
theme(legend.position = "none") +
scale_x_discrete(guide = guide_axis(angle = 45)) +
theme(axis.text.x = element_text(size = 7)) +
scale_fill_tableau()